Abstract: The goal of this study is to look at territory size of the Louisiana Waterthrush during breeding season at three locations; Wall Branch stream located in Rotary Park in Clarksville Tennessee 37040, Dry Creek located in Cheatham Wildlife Management Area in Ashland City Tennessee 37015, Big Hollow located in Beaman Park in Ashland City Tennessee 37015. Mapping these territories in these streams will show not only the territories of the Louisiana Waterthrush present but also possible areas of concentrated use that may aid in the location of their nesting sites.

We will be making KDE maps so we will need to register with google maps to use their satelite version of our map area. The link below is where a lot of the following code came from. https://cfss.uchicago.edu/notes/raster-maps-with-ggmap/ You must register for google maps to get an AI (I used the one year free trial for this and they don’t auto charge once the year is over). OpenMapSource is another option as well if your java supports this. Once you have registered create your base map using google maps address or coordinates.

register_google("AIzaSyDKCi-txHXwSLjB-UURgX1eLzhNVcsXMRg")
basemap<- get_map(location = "1539 Dry Creek Rd
Ashland City, TN 37015", zoom=15)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=1539%20Dry%20Creek%20Rd%0AAshland%20City,%20TN%2037015&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1539+Dry+Creek+Rd%0AAshland+City,+TN+37015&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
ggmap(basemap)

Next, we can put our LOWA encounter cooridinates on the basemap we just created.

DCmap<-ggmap(basemap)+ geom_point(data=DryCreek, aes(Longitude, Latitude)) +
  geom_point(data = DryCreek, aes(Longitude, Latitude), size = 0.1) +
  theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_bw()
plot(DCmap)
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

Now it’s time to create KDE and MCP, but first we have to format a few things. It will be important to first put your data into SpatialPoints we will be using CRS (coordinate reference system) which is the format for most GPS data. The link below is where a lot of the following code came from. https://mhallwor.github.io/_pages/activities_GenerateTerritories 

DCpts <- sp::SpatialPoints(coords = cbind(DryCreek$Longitude, DryCreek$Latitude),
                                proj4string = sp::CRS("+init=epsg:4326"))
head(DCpts)
## class       : SpatialPoints 
## features    : 1 
## extent      : -87.18876, -87.18876, 36.18876, 36.18876  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## put all the data into a single SpatialPointsDataFrame (spdf)
DC_spdf <- sp::SpatialPointsDataFrame(DCpts, DryCreek)
head(DC_spdf)
## # A tibble: 6 x 4
##   Name      Longitude Latitude Number
##   <chr>         <dbl>    <dbl>  <dbl>
## 1 Dry Creek     -87.2     36.2      1
## 2 Dry Creek     -87.1     36.2      2
## 3 Dry Creek     -87.1     36.2      3
## 4 Dry Creek     -87.1     36.2      4
## 5 Dry Creek     -87.1     36.2      5
## 6 Dry Creek     -87.1     36.2      6
#although we don't have seperate species we need to let R know that all points are for one stream system and this will create a list of one
DC_sep <- split(x = DC_spdf, f = DC_spdf$Name, drop = FALSE)

Time to make the MCP.

DCmcp <- lapply(DC_sep, FUN = function(x){rgeos::gConvexHull(x)})

##this makes polygon from the list of one when we told R that all variables are from one stream system
DCmcp <- mapply(DCmcp, names(DCmcp), 
                  SIMPLIFY = FALSE,
                  FUN = function(x,y){x@polygons[[1]]@ID <- y
                  return(x)})

DCmcp <- do.call(rbind,DCmcp)
DCmcp <- SpatialPolygonsDataFrame(Sr = DCmcp,
                                   data = data.frame(Bird = names(DCmcp)),
                                   match.ID = FALSE)
plot(DCmcp)

Now that we have an MCP for the stream, onto KDE.

## Step one: do least squares cross-validation to estimate bandwidth (you may get a warning message but keep going)
bw <- lapply(DC_sep, FUN = function(x){ks::Hlscv(x@coords)})
## Warning in ks::Hlscv(x@coords): Data contain duplicated values: LSCV is not
## well-behaved in this case
## Step two: generate kde

DC_kde <-mapply(DC_sep,bw,
                     SIMPLIFY = FALSE,
                     FUN = function(x,y){
                       raster(kde(x@coords,h=y))})
# This code makes a custom function called getContour. 
# Inputs:
#    kde = kernel density estimate
#    prob = probabily - default is 0.95

getContour <- function(kde, prob = 0.95){
  # set all values 0 to NA
  kde[kde == 0]<-NA
  # create a vector of raster values
  kde_values <- raster::getValues(kde)
  # sort values 
  sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
  # find cumulative sum up to ith location
  sums <- cumsum(as.numeric(sortedValues))
  # binary response is value in the probabily zone or not
  p <- sum(sums <= prob * sums[length(sums)])
  # Set values in raster to 1 or 0
  kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
  # return new kde
  return(kdeprob)}

DC_95kde <- lapply(DC_kde,
                    FUN = getContour,prob = 0.95)

These next plots put MCP on top of KDE 

plot(DC_kde[[1]])+
plot(DCmcp[1,],add = TRUE)

## integer(0)

Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html 

KDEmap<-DCmap +
  stat_density2d(aes(x = Longitude, y = Latitude, fill = ..level..,alpha=..level..), bins = 10, geom = "polygon", data = DryCreek) +
  scale_fill_gradient(low = "red", high = "green")+
  ggtitle("Dry Creek, TN")
Compleated KDE Heatmap

Compleated KDE Heatmap

Big Hollow, TN

Big_Hollow <- read_csv("/Users/Kirstenglish/Desktop/Big Hollow.csv")
## Parsed with column specification:
## cols(
##   Name = col_character(),
##   Latitude = col_double(),
##   Longitude = col_double(),
##   Number = col_double()
## )
View(Big_Hollow)

Just like before, retrieve a basemap of your site location.

basemap<- get_map(location = "Ridgetop Trail
Ashland City, TN 37015", zoom=15)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=Ridgetop%20Trail%0AAshland%20City,%20TN%2037015&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Ridgetop+Trail%0AAshland+City,+TN+37015&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
ggmap(basemap)

Put our LOWA encounter cooridinates on the basemap we just created.

#time to put points on the map
BHmap<-ggmap(basemap)+ geom_point(data=Big_Hollow, aes(Longitude, Latitude)) +
  geom_point(data = Big_Hollow, aes(Longitude, Latitude), size = 0.1) +
  theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_bw()
plot(BHmap)

Put data into spatial points like before to be able to create the KDE and MCP. we will be using CRS (coordinate reference system) which is the format for most GPS data.

Time to make the MCP.

BHmcp <- lapply(BH_sep, FUN = function(x){rgeos::gConvexHull(x)})

##this makes polygon from the list of one when we told R that all variables are from one stream system
BHmcp <- mapply(BHmcp, names(BHmcp), 
                SIMPLIFY = FALSE,
                FUN = function(x,y){x@polygons[[1]]@ID <- y
                return(x)})

BHmcp <- do.call(rbind,BHmcp)
BHmcp <- SpatialPolygonsDataFrame(Sr = BHmcp,
                                  data = data.frame(Bird = names(BHmcp)),
                                  match.ID = FALSE)
plot(BHmcp)

Now that we have an MCP for the stream, onto KDE.

bw <- lapply(BH_sep, FUN = function(x){ks::Hlscv(x@coords)})
## Step two: generate kde

BH_kde <-mapply(BH_sep,bw,
                SIMPLIFY = FALSE,
                FUN = function(x,y){
                  raster(kde(x@coords,h=y))})

# This code makes a custom function called getContour. 
# Inputs:
#    kde = kernel density estimate
#    prob = probabily - default is 0.95

getContour <- function(kde, prob = 0.95){
  # set all values 0 to NA
  kde[kde == 0]<-NA
  # create a vector of raster values
  kde_values <- raster::getValues(kde)
  # sort values 
  sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
  # find cumulative sum up to ith location
  sums <- cumsum(as.numeric(sortedValues))
  # binary response is value in the probabily zone or not
  p <- sum(sums <= prob * sums[length(sums)])
  # Set values in raster to 1 or 0
  kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
  # return new kde
  return(kdeprob)}

BH_95kde <- lapply(BH_kde,
                   FUN = getContour,prob = 0.95)

These next plots put MCP on top of KDE 

plot(BH_kde[[1]])+
  plot(BHmcp[1,],add = TRUE)

## integer(0)

Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html

Compleated KDE Heatmap

Compleated KDE Heatmap

Wall Branch, TN

Create the basemap for the site location. This site is Wall Branch stream in Tennesse

basemap<- get_map((location = " 2561 Alex Overlook Way
Clarksville, TN 37043, 36.497818, -87.267428"), zoom=15, maptype = "satellite")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=%202561%20Alex%20Overlook%20Way%0AClarksville,%20TN%2037043,%2036.497818,%20-87.267428&zoom=15&size=640x640&scale=2&maptype=satellite&language=en-EN&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=2561+Alex+Overlook+Way%0AClarksville,+TN+37043,+36.497818,+-87.267428&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
ggmap(basemap)

LOWA encounter cooridinates on the basemap we just created.

#time to put points on the map
WBmap<-ggmap(basemap)+ geom_point(data=Wall_Branch, aes(Longitude, Latitude)) +
  geom_point(data = Wall_Branch, aes(Longitude, Latitude), color= "black", size = 0.1) +
  theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_classic()
plot(WBmap)
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

WBpts <- sp::SpatialPoints(coords = cbind(Wall_Branch$Longitude, Wall_Branch$Latitude),
                           proj4string = sp::CRS("+init=epsg:4326"))
head(WBpts)
## class       : SpatialPoints 
## features    : 1 
## extent      : -87.26923, -87.26923, 36.49896, 36.49896  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
WB_spdf <- sp::SpatialPointsDataFrame(WBpts, Wall_Branch)
head(WB_spdf)
## # A tibble: 6 x 4
##   Name        Latitude Longitude Number
##   <chr>          <dbl>     <dbl>  <dbl>
## 1 Wall Branch     36.5     -87.3      1
## 2 Wall Branch     36.5     -87.3      2
## 3 Wall Branch     36.5     -87.3      3
## 4 Wall Branch     36.5     -87.3      4
## 5 Wall Branch     36.5     -87.3      5
## 6 Wall Branch     36.5     -87.3      6
WB_sep <- split(x = WB_spdf, f = WB_spdf$Name, drop = FALSE)

Time to make the MCP.

WBmcp <- lapply(WB_sep, FUN = function(x){rgeos::gConvexHull(x)})

##this makes polygon from the list of one when we told R that all variables are from one stream system
WBmcp <- mapply(WBmcp, names(WBmcp), 
                SIMPLIFY = FALSE,
                FUN = function(x,y){x@polygons[[1]]@ID <- y
                return(x)})
WBmcp <- do.call(rbind,WBmcp)
WBmcp <- SpatialPolygonsDataFrame(Sr = WBmcp,
                                  data = data.frame(Bird = names(WBmcp)),
                                  match.ID = FALSE)
plot(WBmcp)

Now that we have an MCP for the stream, onto KDE.

bw <- lapply(WB_sep, FUN = function(x){ks::Hlscv(x@coords)})

## Step two: generate kde
WB_kde <-mapply(WB_sep,bw,
                SIMPLIFY = FALSE,
                FUN = function(x,y){
                  raster(kde(x@coords,h=y))})
# This code makes a custom function called getContour. 
# Inputs:
#    kde = kernel density estimate
#    prob = probabily - default is 0.95

getContour <- function(kde, prob = 0.95){
  # set all values 0 to NA
  kde[kde == 0]<-NA
  # create a vector of raster values
  kde_values <- raster::getValues(kde)
  # sort values 
  sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
  # find cumulative sum up to ith location
  sums <- cumsum(as.numeric(sortedValues))
  # binary response is value in the probabily zone or not
  p <- sum(sums <= prob * sums[length(sums)])
  # Set values in raster to 1 or 0
  kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
  # return new kde
  return(kdeprob)}

WB_95kde <- lapply(WB_kde,
                   FUN = getContour,prob = 0.95)

These next plots put MCP on top of KDE 

## These next plots put MCP on top of KDE, (territory for this map is very linear)
plot(WB_kde[[1]])+
  plot(WBmcp[1,],add = TRUE)

## integer(0)

Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html 

WBkde<-WBmap +
  stat_density2d(aes(x = Longitude, y = Latitude, fill = ..level..,alpha=..level..), bins = 10, geom = "polygon", data = Wall_Branch) +
  scale_fill_gradient(low = "red", high = "green")+
  ggtitle("Wall Branch, TN")
Compleated KDE Heatmap

Compleated KDE Heatmap

Species Distrubution Map Louisiana Waterthrush

library(dismo)
## 
## Attaching package: 'dismo'
## The following object is masked from 'package:ggmap':
## 
##     geocode
library(rgbif)
library(rdryad)
## Registered S3 method overwritten by 'crul':
##   method                 from
##   as.character.form_file httr
library(readxl)
library(ggridges)
## 
## Attaching package: 'ggridges'
## The following object is masked from 'package:ggplot2':
## 
##     scale_discrete_manual
library(viridis)
## Loading required package: viridisLite
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
extent <- extent(-102,-80,10,45)

LOWA_dismo <- gbif("seiurus", species = "motacilla", ext = extent,
             geo = TRUE, sp = TRUE, download = TRUE,
             removeZeros = TRUE)
## 938 records found
## 0-300-600-900-938 records downloaded
LOWA_xy <- as.data.frame(cbind(LOWA_dismo@coords[,1],LOWA_dismo@coords[,2]))
colnames(LOWA_xy) <- c("longitude","latitude")

us <- map_data("world")
ggplot(data = LOWA_xy, aes(x=longitude, y=latitude)) +
  geom_polygon(data = us, aes(x=long, y = lat, group = group),
               fill = "white", color="black") +
  geom_point() + xlab("Longitude") + ylab("Latitude") +
  coord_fixed(xlim = c(-102,-80), ylim = c(10,45))

bioclim <- getData(name = "worldclim", res = 2.5, var = "bio")

names(bioclim) <- c("Ann Mean Temp","Mean Diurnal Range","Isothermality","Temperature Seasonality","Max Temp Warmest Mo","Min Temp Coldest Mo","Ann Temp Range","Mean Temp Wettest Qtr","Mean Temp Driest Qtr","Mean Temp Warmest Qtr","Mean Temp Coldest Qtr","Annual Precip","Precip Wettest Mo","Precip Driest Mo","Precip Seasonality","Precip Wettest Qtr","Precip Driest Qtr","Precip Warmest Qtr","Precip Coldest Qtr")

bio_extent <- extent(x = c(
  min(LOWA_xy$longitude),
  max(LOWA_xy$longitude),
  min(LOWA_xy$latitude),
  max(LOWA_xy$latitude)))

bioclim_extent <- crop(x = bioclim, y = bio_extent)
bioclim_model <- bioclim(x = bioclim_extent, p = LOWA_xy)
presence_model <- dismo::predict(object = bioclim_model, 
                                 x = bioclim_extent, 
                                 ext = bio_extent)
gplot(presence_model) + 
  geom_raster(aes(fill=value)) +
  geom_polygon(data = us, aes(x= long, y = lat, group = group),
               fill = NA, color="black") +
  geom_point(data = LOWA_xy, aes(x = longitude, y = latitude), color = "black", alpha = 0.5) +
  scale_fill_gradientn(colours=c("brown","yellow","darkgreen"), "Probability") +
  coord_fixed(xlim = c(-102,-80), ylim = c(10,45)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Probability of LOWA Occurrence") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + theme(legend.position = "right") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())